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Variational methods with coupled Gaussian functions for Bose-Einstein condensates 
with long-range interactions. I. General concept 

Stefan Rau, Jorg Main, and Giinter Wunner 
InsUtut fur Theoretische Physik 1, Universitdt Stuttgart, 70550 Stuttgart, Germany 

(Dated: August 12, 2010) 

The variational method of coupled Gaussian functions is applied to Bose-Einstein condensates 
with long-range interactions. The time-dependence of the condensate is described by dynamical 
equations for the variational parameters. We present the method and analytically derive the dy- 
namical equations from the time- dependent Gross-Pitaevskii equation. The stability of the solutions 
is investigated using methods of nonlinear dynamics. The concept presented in this paper will be 
applied to Bose-Einstein condensates with monopolar 1/r and dipolar 1/r interaction in the sub- 
sequent paper [S. Rau et al., Phys. Rev. A, submitted], where we will present a wealth of new 
phenomena obtained by using the ansatz with coupled Gaussian functions. 

PACS numbers: 67.85.-d, 03.75.Hh, 05.30.Jp, 05.45.-a 
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I. INTRODUCTION 

The experimental realization of Bose-Einstein conden- 
sates (BEG) with ^^Cr atoms [1, 2], with a strong dipole- 
dipole interaction has given new impetus to theoretical 
investigations of BEG with long-range interactions. 

The theoretical description of Bose-Einstein conden- 
sates in the dilute limit in the framework of the ex- 
tended Gross-Pitaevskii equation (GPE) is well known. 
The derivation of the extended Gross-Pitaevskii equa- 
tion from a many particle Schrodinger equation is part of 
many text books on quantum mechanics or Bose-Einstein 
condensates [3] . For a long-range interaction of the form 
W\r{r,r') (X \r — r'\", the time-dependent GPE can be 
brought into particle number scaled dimensionless form 
using appropriate units (for monopolar or dipolar con- 
densates, see [4, 5]), and reads 



- A + jIx^ + ^y + 7f z^ + Sttgsc mr, t)f 



+ /dVWu.(r,r')|V(r',i)r 



^(r) = i-V'(r,i). (1) 



The terms in Eq. (1) describe the short-range contact 
interaction between two particles, the s-wave scattering, 
14c = 87rasc|'0('")P: 8- harmonic model for external mag- 
netic trapping of the condensate, Vt = 7^x^-1-7^2/^-1-7^2^, 
and long-range interactions between two particles 

Vi,{r) = f d^r'Wi,{r,r')\ij{r')\^ . 



The mean field energy reads 

E^f - (tAI - A + K + (Vsc + Mr)/2|V^) 



(2) 



So far, in most publications one of two methods for solv- 
ing the GPE is used. The ground state for the GPE with 
long-range interaction in Eq. (1) is obtained by minimiz- 
ing the energy functional (2) using different approaches. 



The first method are numerical lattice calculations [6] , 
either the minimization of the energy with conjugate gra- 
dients or imaginary time evolution of an initial wave func- 
tion using the split-operator method and EFT. The nu- 
merical calculations are on the one hand very accurate 
if they are carried out on sufHciently large grids, on the 
other hand, however, they may turn out laborious and 
may take a long computational time. 

The second well-established method is using a sim- 
ple variational ansatz. A common ansatz is to assume 
a Gaussian type wave function [7-11]. This technique al- 
lows to gain physical insight, as it often provides qualita- 
tively correct, although quantitatively inaccurate results. 
The extension and improvement of variational techniques 
for BECs with long-range interactions is the major chal- 
lenge of this paper. 

We propose, as a third approach, an improved vari- 
ational ansatz with coupled Gaussian functions. The 
method was originally proposed by Heller [12, 13] to de- 
scribe atomic and molecular quantum dynamics. Using 
the ansatz 

tpix, = X! exp{i[(a; - Xt^n)At^n{x - Xt^n) 
n 

+ Pt,n{x - Xt^n) + lt,n]} (3) 

where the symmetric matrix At^n, the vector pt^n and the 
scalar jt,n describe the width, momentum and weight of 
the Gaussian wave packet, respectively. Heller approxi- 
mated the dynamics of quantum wave packets following 
classical trajectories. The method was recently success- 
fully applied to the dynamics of atoms in external fields 
[14, 15] using up to iV = 100 coupled functions. 

We apply this method to Bose-Einstein condensates 
with long-range interactions using an ansatz with N cou- 
pled Gaussian functions centered at the origin, Xt^n — 0, 
and with pt^n ~ 0, viz. 



N 



i^(.r,t)^J2 



fe=l L 



,17 



n (»"■"■■) 



(4) 



where a G {x, y, z} for a BEC without symmetries, a G 
{g, z} for an axisymmetric, and a = r for a spherically 
symmetric BEC. 

The ansatz in Eq. (4) is rather general and is not only 
able to describe condensates with spherical or axial sym- 
metry, but also nonsymmctric condensates in arbitrary 
trap geometries or even anisotropic solitons [16, 17]. The 
method is also well suited for long-range interactions be- 
cause by its very nature it requires integrals over the en- 
tire wave function which show a similar behavior as the 
integrals for local interactions. Note that, as is typical for 
a Gaussian basis set, the individual Gaussian functions 
in Eq. (4) are not orthogonal amongst each other. In this 
paper we present the theoretical concept of the method 
and derive the necessary equations for arbitrary long- 
range interaction. For monopolar (1/r) or dipolar (l/r^) 
interaction selected results have already been presented 
in [18]. Results elaborated in more detail are subject 
of the subsequent paper [19]. With the use of multiple 
Gaussians and thus an extended set of variational param- 
eters, we are not only able to describe the stable ground 
state and the metastable stationary states, but can also 
identify the types of bifurcations where branches emerge 
or stability changes take place. 

The paper is organized as follows. In Sec. II we ap- 
ply a time-dependent variational principle to the Gross- 
Pitaevskii equation (1) and obtain dynamical equations, 
which describe the time-dependence of the variational pa- 
rameters. In Sec. Ill we evaluate the integrals that are 
needed to set up this nonlinear set of dynamical equa- 
tions. In Sec. IV we present different methods for ob- 
taining stationary solutions, and in Sec. V we investigate 
the stability of those states. Gonclusions are given in 
Sec. VI. 



by the respective Schrodinger or Gross-Pitaevskii equa- 
tion 



H^it) 



dt 



tPit) 



(6) 



The variational principle of McLachlan minimizes the dif- 
ference between the left and the right hand side of the 
respective Schrodinger equation (6) with respect to the 
trial wave function. 



mt) - Hi^im^^ 



(7) 



For any t, ipit) is supposed to be fixed and given, and 
the quantity / is minimized by varying (j). Afterwards 
(j) is set equal to <j> — ip. The time dependence of the 
trial wave function carries over to the time dependence of 
the variational parameters, 'ip{t) = ip{z{t)). We consider 
variations of / in Eq. (7) with respect to 0, 

SI = {S(^\(^) + {q)\Sq)) + i{5dp\Hi,) - i{Hyj\S^) 

- iHtp) + {(I) ~ iH^\S4>) , (8) 



where the variation of the time derivative of the wave 
function S(j) carries over to variations of the parameters 
z. 



\s4>)^mz,z)) 



dtp 
dz 



Sz 



dip 
dz 



Sz 



(9) 



The first term vanishes, since we minimize / under the 
condition that '0(i) is fixed and therefore that all param- 
eters z are fixed as well, and we obtain 



M^ 



d f dtp . 
dz \ dz 



dtp_ 
dz 



Sz 



(10) 



II. TIME DEPENDENT VARIATIONAL 
PRINCIPLE 



We insert Eq. (10) in Eq. (8), set (f> — tp, and obtain as 
condition for the vanishing variation of / 



The first variational principle in order to optimize a 
variational wave function that comes to mind is mini- 
mizing the mean field energy functional. For few varia- 
tional parameters the analytical calculation is done eas- 
ily, but it gets increasingly difficult if more variational 
parameters are included in the wave function. For a 
BEC with long-range interaction and more than three 
parameters, it turns out to be almost impossible to cal- 
culate all derivations analytically. Therefore we will in- 
troduce a different approach based on the Dirac-Frenkel- 
McLachlan variational principle [20, 21]. The applica- 
tion of this variational principle yields a set of differential 
equations for the parameters z of the trial wave function 
tp{t) =tp[z{t)), where 



z{t)^{zi{t),Z2{t),...,ZM{t)) 



(5) 



is the vector consisting of all M variational parameters. 
The time dependence of a quantum system is described 



dz 



tp + iHtp ) + {tp- iHtp 



g.i).0.,U, 



The variational parameters z are complex quantities and 
therefore the variations Szk and (5i^ for k = 1, . . . , M are 
independent. Therefore both brackets of Eq. (11) have 
to vanish separately. This finally yields the equation 



dtp 
dz 



itP- HtP) =0 



(12) 



which is easily transformed to an implicit dynamical set 
of equations Kz = —ih for the variational parameters 
with 



K 



dtp 
dz 



dtp 
dz 



h 



dtp_ 
dz 



HtP 



Up to this point there has been no specification of the 



trial wave function -0 or the Hamiltonian H . The varia- 
tional principle can be applied to both, linear and non- 
linear Hamiltonians. In the following we apply the time- 
dependent variational principle (TDVP) to a trial wave 
function given as coupled Gaussian wave functions, and 
derive dynamical equations for the variational parame- 
ters of each Gaussian. 



A. Dynamical equations for condensates without 
symmetries 

We choose a superposition of N Gaussians as trial wave 
function, 



N 



N 



^{r, t)^J29'^J2 ^i^-'+^W+^^'+i") . (13) 



k=l 



k=l 



The a^ for a G {x,2/, z} denote complex Gaussian 
"width" parameters in the three spatial directions, and 
7*^ the complex amplitude/phase parameters. The sys- 
tem is described by the extended GPE (1) with the 
Hamiltonian brought to the scaled "natural" units for 
the respective long-range interaction. 



H = -A + Kfr(r) , 



(14) 



where Vcfi{r) — Vt{r) + T4c(^) + Mr(^) is the sum of 
the trapping, scattering and long-range potential. We 
use Eq. (12) obtained from the TDVP for the complex 
variational parameters 

— ('^1 r^N 1 N 1 N 1 N\ 

Z 1^7 , . . . , 7 ,U^,...,U^,Uy,...,Uy,U^,...,U^), 

(15) 
and first calculate the time derivative '0 of the coupled 
Gaussian wave function. The derivation carries over 
to time derivatives of the Gaussian width and ampli- 
tude/phase parameters. 



N 



N 



Ti^s'^^^' (^'''- + y^^'y + ^'''^ + ^'') s''- 



dt 



(16) 



fc=i k=i 

Second, we apply the Laplace operator of the Hamilto- 



nian in Eq. (14) in Cartesian coordinates, 

N 



+ ay + a^] 



-A0 - ^|-2i K 

k=l *- 

+ 4[(a^)%2 + (4)%2 + (a^)2z2]|/, (17) 

and obtain the complete expression for the kct in Eq. 

(12), 



I — 1 V ^ 



k=l 



jv;ff(r)-2i[a^ + a: 



{al) x^ + ial) v'^ + {a'lfz 



(18) 



The sorting of Eq. (18) according to powers of the coordi- 
nates X, y, z results in a sum of products of a polynomial 
of second order and the Gaussian g*^ , 



N 



i^-Hi,=J2 



'^'o + livL^' 



fe=i 

+ V,'yy' + Vl^z' 



V,s{r) 



with the newly defined quantities 

v^ = -7'= + 2i (a^ + a^ + a';) , 
0^2% = -^ {a-if - ai ; ae{x,y,z} 



(19) 

(20a) 
(20b) 



^ J 



Now, we calculate the derivatives oi ip — J2i=i 9 ^^ 
Eq. (12) with respect to the variational parameters z 
in Eq. (15) for each / = 1, . . . , iV: 






dip 



'^xW ; 



€ {x, y, z} . 



(21) 



Finally, Eq. (12) results in the 4Af-dimensional system of 
equations 



J 



v'g' 



N 

E 



'0 + :, {yi.^^ + ylyv" + yt^^) - y.^(r) 



(22) 



with ^ — 1, . . . , A'^ and ry = 1, x, j/, z, which can be sorted as 

N N N 

E(5' I /K + ^EE(5' I -I I 9')yla = EV I ^eff I /); l = l,...,N; 



fc=l 



fc=l a 



k=l 

N 



N N N 

E(5' I 4 I a'H + 7;T.T.(9' I ^'4 I 9'Mc. - E(3' I 4^off I g'); f3 e {x,y,z}-J ^ 1,. . .,N., 



(23) 



fc=i 



fe=l a 



k=l 



where for a BEC without symmetries a,/3 € {x,y,z}. 
Equation (23) can be written in the form of a matrix 
equation 

Mv = r, (24) 

with the Hcrmitian, positive definite 4iVx 4iV matrix M, 



)ik \^ )ik {y )ik \z )ik 

1\/T - I i^'^)kl {x'^)lk {x'^y'^)lk {x^Z^)lk 

^^ ' {y^)ki iy'x^)ki {y%k {y^z^)iu 

,{z'^)ki {z'^x'^)ki {z'^y'^)ki {z^)ik 



(25) 



where all terms are N x N matrices for fc = 1 , . . . , A^ and 
I = 1, . . . ,N. As an example the term {x'^)ik reads 



{-') 



Ik 



[9'-'\x'\9''=') 



,{g'=''\x^\g'^=') 



l = ll^2l„k=N\ 



\x^\g 



(^gl=N\^2\gk=N) 



(26) 

The terms denoted {y'^)ik,{z'^)ik,{x^y'^)ik, ■ ■ ■, are ana- 
log to the example and obtained by replacing the x^ by 
2/^, z^, x^y^, . . ., respectively. The vectors v and r in 
Eq. (24) are 



V = 



N 

E 

k=l 



2 ^2,x 
2^2,y 



{9'\y'vM9'' 

\{9'\z^V,n\g>^) ) 



where each entry is a vector of length N for k = 1, . . . ,N 
and ^ = 1, . . . , A^ respectively. 

By solving the definitions of (wq , V2) in Eq. (20) for the 
time derivatives of the Gaussian parameters, we obtain 
AN dynamical equations for the Gaussian parameters z, 



2i(a^ 



k\2 



«) 



:V' 



a e {x,y,z}; fc = 1, 



(29a) 

.,N, 
(29b) 



keeping in mind that the quantities (vqjV^) — 
{vq, ¥2^^, V2^y, ¥2^^) constitute the solution vector to the 
linear set of equations (24). These linear equations con- 
tain basic Gaussian integrals in the matrix (25) on the 
left hand side, as well as integrals with the interaction 
terms of the Hamiltonian in the vector (28) on the right 
hand side. The necessary integrals will be calculated an- 
alytically in Sec. Ill for condensates with rather general 
long-range interactions. 



B. Dynamical equations for condensates with axial 
or spherical symmetry 

If the GPE describes a system that is constrained by, 
e.g., axial or spherical symmetry, the results obtained 
in Sec. II A may be adapted to the case at hand. There- 
fore we introduce respective coordinates {g, 4>, z) for axial 
symmetry and (r, 0, 0) for spherical symmetry and choose 
suitable trial wave functions 



N 



^(r,i) = ^ 



fe=i 



=17 



n(^ 



(30) 



where for axially symmetric BEC we have a G {g, z}, and 
for spherically symmetric, e.g., monopolar BEC we have 
a = r. These trial wave functions reduce the number of 
complex parameters to 3A^ (0^,0^,7*^) and 2N (ajr,7'^) 
(k — 1 . . . , N) , respectively, where N is the number of 
Gaussians. 

The procedure is the same as in Sec. II A: First, 
we calculate the time derivative ipit) which can be ob- 
tained by Eq. (16) by simply setting a^ = Uy = a^ and 
a^ — Gy — a'^ — a'^, accordingly. Second, the Laplace op- 
erator is applied to the coupled Gaussian wave function, 
and finally with the respective definitions of the vectors 
{v^i,y2V^^g,y2V^^^)'^ and (wg, 1/2^2';^)'^, the dynamical 



(27) equations can be written as 



2i(2a' 



t = -4(a^)2 - il^^^a; k = l,...,N;ae{g,z}., 



(31) 



(28) for a BEC with axial symmetry {D — 2), and 



7*= = 6ia^ - Vq 



fc = 1, 



,iV, 



(32) 



The quan- 



^ ' k ^ u\rk \ j^g^^g ^p |-,g 



for a BEC with spherical symmetry (JJ — 1 

titles («o^ 1/21^2':^, ^h^lzf and (i;§, 1/2^^2':^) 
calculated from an adapted set of linear equations analog 
to Eq. (23), but for a BEC with axial symmetry, D — 2, 
with 

a,/3e{g,z}, (33) 

and for a BEC with spherical symmetry, D ~ I, with 

a,P = r. (34) 

As in Eq. (24), we can rewrite both respective linear 
sets of equations in the form of a matrix equations 



Mv = r. 



(35) 



with the Hcrmitian, positive definite {D + l)Nx (D+1)N 
matrix M analog to Eq. (25) but with a reduced num- 
ber of blocks, {{g'^)ki,{z^)ki,---} or {(r^)^;, . . .}, and 



the (-D + l)7V-dimensional respective vectors v and r for 
D = 2and D = l. 

We have apphed the time-dependent variational prin- 
ciple of Dirac Frenkel and McLachlan to the extended 
Gross-Pitaevskii equation and a trial wave function with 
coupled Gaussian functions. With the resulting dynami- 
cal equations (29), (31) and (32) in respective symmetries 
we can calculate the time dependence of the wave func- 
tion Tp by calculating the time dependence of the varia- 
tional parameters. To set up the equations (29), (31) and 
(32) we have to solve the set of linear equations for the 
quantities (wq , V2,a) , k — 1, . . . ,N;a — x,y,z {a = g,z 
and a — r respectively) in Eq. (23). All integrals of the 
matrix and the right hand side are calculated analyti- 
cally. 



III. CALCULATION OF THE INTEGRALS 

For clarity we sort all integrals as they appear in this 
matrix equation: the integrals for the matrix (25) in 
Sec. Ill A, and the integrals of the right hand side (28) 
in Sec. IIIB. Then we calculate the mean field energy 
and the chemical potential for a BEG with long-range 
interaction in Sec. IIIC. 



A. Computation of the matrix M 

The integrals of the matrix M in Eq. (25) are all 
of the form {g^\g''), {g^\xa\g''), and (g'|a;^a;||g'^), with 
Xa,x^ G {x, y, z}, and each Gaussian function g'^ defined 
in Eq. (13). All integrals are easily calculated from the 
simplest integral {g^\g^) with the use of the relation, 



(g'|xf4'^l/|/) = (-i) 



\+u 



d' 



d" 



d{a^rd{a''p) 



A9'\V\g'), 
(36) 



with A, i^ = 0, 1, 2, . . . , and V an arbitrary potential. 

To facilitate reading and to shorten the extensive terms 
in integrals or integral solutions, we introduce the follow- 
ing abbreviations: 



C=aa+«L', ae{x,y,z}, 
yj ==7*- (7^)% 



(37a) 
(37b) 
(37c) 
(37d) 
(37e) 
(37f) 



We start with 

00 
(g'|/)=e'^" j 0^'^'dx j c<y"dy 

3/2 i'v'"' 



V-i«J;'\/-i<V" 



-,kl 



(38) 



where, as can be seen easily from the definition of the 
Gaussian trial wave function (4), the imaginary parts 
of the widths are to remain positive. Therefore the 
imaginary parts of the occurring combinations also fulfill 
Im [a'= - a'*] > 0. 

The application of the relation (36) to the norm inte- 
gral (38) provides with A = 1, j^ = the integrals 



(5' I ^' I /) 



{9' I y 



'1^2 



(3' 

and with A 



(5' I ^' 



(5' I y 



{9' I z' 



and 

(5' I ^V 



{9 



' I x^z^ 



J I ..2^2 



/I 2 

\9 I y 



V^pii" 



') = 



2(-i4o'/'v^y=i^ 



2V^i^(-i<)''V^ 



/) 



V 



'/" / • ki 



4 (-ia^O'^' y^-ia^V-ia^ 



2v/^kFy^(-ia^0'/^' 
- V = 2 the integrals 

9') = 



(39a) 

(39b) 
(39c) 



37r-V^e'^" 
4y^i^y^(-iaJ0''/^' 



4(-i40'^^(-i<)'^'V^i^' 



(40a) 

(40b) 
(40c) 

(41a) 



V 



TT '"e 






, (41b) 



^2, 



4y=kF(-ia^0'' (-ia^O 



.3/2- 



(41c) 



Since the matrix is now complete, we turn to the more 
challenging integrals needed for the right hand side of 
Eq. (24), i.e., the vector r defined in Eq. (28). 



B. Computation of the vector r 

The right hand side of Eq. (24) contains the trapping 
term, Vt, the scattering term, Vsc, as well as the more 



complicated long-range interaction term, Vir. In the fol- 
lowing equations the g^ defined in Eq. (13) represent the 
individual Gaussians of the trial wave function. 



1 . Trapping 
We start with the term of the trapping potential 

TrV^e"'" (7^aJ'a^' +7^a«af +7fa^'a«) 



(42) 



For the right hand side of Eq. (24) we also need the terms 

with Xq, (E {x,j/, 2;}. They are directly obtained with the 
help of Eqs. (40) and (41). 

2. Scattering 

The second interaction term of Vcs on the right hand 
side of Eq. (24) contains the nonlinear contact interac- 
tion Vsc — 87rasc |'0('")l of the s-wave scattering. Follow- 
ing the same procedure as above, we start by calculat- 
ing (g'JKiclff'^), before obtaining the terms (5'|a;^T4c|ff'^) 
with Xa G {x, y, z} using the relation (36). 

(.9' I ■l^c I 9") - (.9' I StTOsc 1^1' I /) 
N ,. 

== ^ / dV (8^a.e 9\r)* g' {r)* g\r)g\r) 



N 



SascTi-"/' ^ 



. ., / '• klij / • klij / '• klij 

,_;^i y ~\ax ' y -lay ' y -\az 



(43) 



as long as Ima^'*^ > for a e {x, y, z}, which is true 
[see Eq. (37) and note that Ima^ > for all width pa- 
rameters]. Again we use relation (36) and get the three 
remaining integrals, 

{9' I xV.e I /) 



*j'^i —\a 



klij 



3/2 



klij 



I • klij 



:, (44a) 



{9' I y'v^c I /) 



N 



- 4a.c^''/^ Y. 



ni7 



- , / • klij ( ■ kiij\''^ r 

(3' I ^'^sc I /) 



klij 



-, (44b) 



JV 



^ 4a.c^''/^ ^ 



Ai7'- 



- . T / ^ HIT 
2,j-i y—iux 



klij 



^kuA 



3/2- 



(44c) 



3. Long-range interaction 

The most challenging calculation surely is that of 
the long-range interaction term (see Eq. (1)) V\i ~ 
/d3r'Wir(r-r')|^(r')|', 

{9' I Vir I /) = E Z'^'''/^'''' ^'^^'^ " ""'^ 



«,i=i 



X.9-' (r')3*(0 5' W/(r). (45) 

Introducing relative (rel) and "center of mass" coordi- 
nates (cm) via 



1 / 1 A frrcl 

2 l-l 1/ Ucm 



(46) 



and keeping the Jacobian determinant (i/s) in mind, 
Eq. (45) is transformed to 



N -.^kiij 



{9' I Mr I /) = ^ ^^ Jd^r,,i Jd^r.^WiArrei) 

X exp< - f a^' (Xcm + aJrol)^ + Oy' (ycm + Viclf 

+ a^-' (j/cm - Z/rol)^ + a'/ (Zcm - Zrcl)^ j ^ , (47) 

with the abbreviations (37). We use the integral 



j{a{x+uf+b{x-uf)^^ 



v/-i(a + fo) 



e'^TT (48) 



which, for Im[a + b] > 0, is easily transformed into a 
standard Gaussian integral after completing the square 
in the exponent, and solve the center of mass integral in 
Eq. (47) 



TV 



(g' I Mr I /) = J2 -^ 



'/2. 



-,17 



klij klij klij 

lax ay ■'uz 



-rklij 
^rel ' 



(49) 



with 



rklij 



dV. 



Wlr(T-rcl) 



X exp < 1 



-^ij„kl 
^x "ic 2 



ax 



klij 



^rcl 



rklij 



"■y % „2 

S/rcl 



ij „kl 



ajfa 



a 



klij 
V 



z "z 2 



(50) 



The relative integral Ij.^[ depends on the specific form of 
the two-particle interatomic interaction iyir('''rei)- The 
integral obtained with the coupled Gaussian method is 



formally the same as for the calculation with a single 
Gaussian trial wave function, except for the complex co- 
efficients in the exponent, a'^Ja^/a^^^ for a e {x,y,z}. 
Therefore the method of coupled Gaussian functions 
is applicable to all two-particle long-range interactions 
which can be solved with the simple single Gaussian wave 
function. 

As an example, we present the results for the relative 
integral for dipolar interaction 



Wlr(r-rcl) = Wd(rrcl) 



1 



kvd^ 



r-roi 



(51) 



where I^^i can be expressed in terms of elliptic integrals 

[16, 17], 

4:7T 

Irol = -Z- [KxHyRo {^1,1^1, l) - l] • (52) 

Kx, Hy are complex combinations of the Gaussian widths. 



axa'raz 



klij ij uj 

ay ■'azaf 

ij f^i klij ' 

"v " 



\ ay^K-Uz 



and Rd is the elliptic integral of the second kind in Carl- 
son form [22, 23], 



RD{x,y,z) 



dt 



27 ^{x + t){y + t){z + tf 



Numerically it is convenient to use Carlson's formulation 
for elliptic integrals because there are very fast converg- 
ing algorithms available [22, 23] even for complex argu- 
ments x, J/, z G C 

The three additional integrals (g'lx^Mrjs'^) for Xa <= 
{x,y,z} needed to complete Eq. (24) are obtained us- 
ing derivatives of (g'j Mrj'?'^) with respect to the Gaussian 
widths a^, see Eq. (36), 

{9' I -iv. I /) = -1^(5' I V, I /) . 



calculate the kinetic term and apply the Laplace operator 
to the coupled Gaussian wave function. 



(5'|A|/ 



-27rV^e'^'° 
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W,kl 



(54) 



Now, all integrals that are needed for the linear equa- 
tions for the quantities {vqtV^o), k = l,...,iV;a = 
X, y, z in Eq. (23), as well as all integrals for the mean field 
energy and the chemical potential have been calculated. 
We are able to set up the dynamical equations for the 
Gaussian parameters. These dynamical equations can 
be solved using three different methods, either by mini- 
mization of the mean field energy, by the search for fixed 
points of the dynamical equations, or by imaginary time 
evolution of an initial wave function. 



IV. COMPUTATION OF THE GROUND STATE 
AND EXCITED STATES 

There are three different methods available to calculate 
variational solutions. The solutions obtained via mini- 
mization of the mean field energy or via evolution of an 
initial wave function in imaginary time, are limited to 
the stable ground state. The method of a highly non- 
linear root search of the dynamical equations (29) yields 
all stationary states of the GPE, the ground state and 
collectively excited states. 

The sensitivity of the methods on the initial values 
greatly differs. While the minimization and the imagi- 
nary time evolution are relatively robust, the root search 
requires sufficiently accurate parameters, especially for 
an increasing number of variational parameters. For 
more coupled Gaussian functions, however, the mini- 
mization of the mean field energy and the imaginary 
time evolution get increasingly time-consuming. Results 
should therefore be obtained with the nonlinear root 
search, the other two routines should only be used to 
calculate appropriate initial values. 



C. Energy functional and chemical potential 



A. Minimization of the mean field energy 



We calculate the mean field energy and the chemical 
potential, 

N 1 



k,l=l 

N 

/i= ^(y'l-A + yt + V^c + Mrl/ 



(53b) 



The terms for trapping, scattering, and the long-range 
interaction have been evaluated in Sec. IIIB. We now 



One method of obtaining the ground state of the Gross- 
Pitaevskii equation is to minimize the mean field energy 
functional in Eq. (2). For multiple coupled Gaussians, 
the analytical calculation of all derivatives with respect 
to the variational parameters, e.g., for calculating the 
gradient, is not possible. Therefore we use a numerical 
minimization routine that uses the energy function values 
only. For an increased number of Gaussian parameters, 
the accuracy of this method is limited, but the results 
may be used as initial values for a root search of the 
dynamical equations, which provides much more reliable 
and accurate results. 



B. Root search for fixed points of tlie dynamical 
equations 

The full three-dimensional calculation for condensates 
with long-range interaction includes 4iV complex vari- 
ational parameters of the wave function. We search for 
these solutions of the nonlinear dynamical equations (29), 
(31), and (32) where the dynamic is trivial, i.e., the sta- 
tionary states. The phase of all Gaussians is defined by 
the chemical potential (see Eq. (55)). 

We use a wave function with N complex amplitude and 
phase parameters 7*^ , . . . , 7^ . Since the wave function 
has to be normalized, we have to ensure {ip\4') = 1- Iii 
summary, wc have to find roots of the following system 
of equations 

a^ = -4(a^)'-^^2';.-0, 



4(4)' -0^2^-0, 



-4(a^) 
2i(«^ 



2 
2^2,-""' 



-M. 



(V'|^)-i = o, 



(55) 



for k = 1, . . . ,N. The quantities vo and V2,x,y,z consti- 
tute the solution vector to the set of equations (23). The 
right hand side (28) of this set of equations contains the 
integrals of contact, trap and long-range interaction. The 
calculation of the dipolar interaction yields elliptic inte- 
grals which are evaluated with the help of fast converging 
algorithms [22, 23]. 

The root search itself is highly nonlinear and may be 
performed with an algorithm based on the Powell hy- 
brid method [24] . The success and the speed of any nu- 
merical root search greatly depends on the number of 
variational parameters. Therefore we reduce the num- 
ber of parameters prior to this routine as much as pos- 
sible. Since the ground state of e.g., dipolar conden- 
sates in an axially symmetric trap is axially symmetric, 
a^ — Qy = a'^, it is possible to reduce the first 3A^ equa- 
tions in Eqs. (55) to 2N equations. For condensates with 
l/r interaction, the wave function as well as the station- 
ary states are spherically symmetric. Therefore, we are 
able to further reduce the number of width parameters 
using a^ = Uy = a'^ = a'^. 



C. Imaginary time evolution 

The third method for finding solutions of the GPE is 
the imaginary time evolution of an initial wave function. 
Although the calculation for a distinct scattering length 
may take a long time, the routine is rather insensitive 
to the choice of the initial parameter values, even for a 
large number of variational parameters. Therefore, the 
imaginary time evolution is very useful in the context of 



the calculation of very accurate input values for the root 
search. 

We substitute t ^ r = it in the GPE and calculate 
the (imaginary) time evolution of the dynamical equa- 
tions. For the linear Schrodingcr equation this leads to a 
damping of all states with a factor according to their re- 
spective energy eigenvalues cxp(— i^nr). Sufficiently long 
integration with respect to the imaginary time coordinate 
and renormalizing yields the ground state of the Hamil- 
tonian. The method can also be applied to the nonlinear 
Gross-Pitaevskii equation. 

In a next step, we linearize the dynamical equations 
in the vicinity of the fixed points in order to analyze the 
stability and possible bifurcations. We can also calculate 
fluctuations dip of the stationary wave function. 



V. STABILITY PROPERTIES OF THE 
VARIATIONAL FIXED POINTS 

The standard method for analyzing the stability of so- 
lutions of the GPE is to perturb the wave function of 
the solution and linearize the GPE, which leads to the 
Bogoliubov-de Gennes equations [3]. In this section, we 
make a different approach and apply methods of nonlin- 
ear dynamics to analyze the stability of the condensates. 

The application of the TDVP to a wave function con- 
sisting of coupled Gaussian functions in Sec. II led to a set 
of dynamical equations. In Sec. IV we described meth- 
ods for searching for stationary solutions of the equations 
(55). Minimization of the mean field energy, imaginary 
time evolution or root search yield variational parame- 
ters z^^ for stationary solutions of the GPE, so called 
"fixed points" (FP). In the case of N coupled Gaussian 
functions these parameters are 3A'^ complex widths, TV for 
each a^,a^,a'^ and N complex amplitudes/phases 7*^: 



vFP 



(7 






FP. 



1, 



,N. 



(56) 



In the case of an axially or spherically symmetric BEC, 
the procedure that follows is the same, but with reduced 
sets of subscripts, {g, z) or even (r) for the width param- 
eters. 

We will investigate the stability of these stationary so- 
lutions with the help of a perturbation of the parameters 
at the fixed point. If the parameter set for the stationary 
solution is indeed the minimum of the mean field energy, 
we expect the solution to be stable. In this case small 
changes of the variational parameters will only result in 
a quasi-periodic oscillation confined to the vicinity of the 
original fixed point. By contrast, as we will see, if the 
stationary fixed point is a hyperbolic fixed point, small 
perturbations will lead to an exponential growth of the 
solution. For the following calculation of the stability 
matrix, it is irrelevant which method was used to ob- 
tain the parameters of the stationary solution in the first 
place. 

To observe the time dependence of the perturbations, 
we use the dynamical equations for the Gaussian param- 



eters (29). The fixed point obviously fulfills Eqs. (55), 

7'=(^"^)--M, 

d^(^FP)=0, (57) 

for a — X, y, z; fc — 1, . . . , A^. For small deviations from 
the fixed point, we first split the equations above into 
real (Re) and imaginary (Im) parts (indicated with the 
tilde z — > 5) in order to linearize them, 



6z = JSz. 



(58) 



6z denotes the deviation of the variational parameters 
from those at the fixed point, z = z^^ + Sz, and J 
denotes the 8A^ x 8A^-dimensional real valued Jacobian 
matrix at the fixed point 



J = 



f) {.,k,Kc -.fc.Ini ^k.Kc ^k,lYn\ 



(59) 



with a, P — x,y, z and k,l — 1, . . . , N. The eigenvalues 
of J determine the characteristic stability of the fixed 
point in whose surroundings the linearization takes place. 
In the coordinates given by the eigenvectors of J, all 
differential equations take the form 



which have the simple solution 



,8iV, 



(60) 



(61) 



The eigenvalues occur in pairs, i.e., if Xi is an eigenvalue, 
—Xi is also an eigenvalue of J. Since the Jacobian matrix 
is not symmetric, the eigenvalues A, are real or complex, 
and there are two possibilities. If all eigenvalues Xi are 
purely imaginary the time evolution of the perturbation 
remains confined, since the solution (61) is oscillating. In 
contrast, if at least one of the real parts of the eigenvalues 
(ReAi) is non-zero, any perturbation in the direction of 
the corresponding eigenvector will grow exponentially. 

Using this method we observe the behavior under small 
perturbations. We are able to investigate the stability of 
any fixed point that we obtain for example from a root 
search of the dynamical equations. 

For unstable fixed points the methods of nonlinear dy- 
namics also allow us to gain insight in the collapse mech- 
anism by analyzing variations of the wave function char- 
acterized by the respective eigenvectors 



r- fc.Rc c k,lm r fc,Rc 



Sa 



fc,In 



(62) 



with k — l,...,N; a G {x,y,z}. If this eigenvector 
6zi is axially symmetric, i.e., Sa^^ = Say^ for all k, the 
perturbation of the condensate is symmetric. If, how- 
ever, the eigenvector breaks the axial symmetry, i.e., 
Sa^ j = —Stty j for all k, the perturbation leads to an 



asymmetric oscillation or collapse of the condensate, de- 
pending on whether the respective eigenvalue is purely 
real or imaginary. 

With the variations of the variational parameters of 
the eigenvector we calculate the expansion of the wave 
function at the fixed point and get the solution for the 
linearized variation of the wave function 



SMr,t) 
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Ait 



(63) 



with a e {x,y,z}. The respective eigenvalue is denoted 
Xi and g'^l^^ is the unperturbed Gaussian k at the fixed 
point. 

Is the approach based on the stability eigenvalues 
of the Jacobian fully equivalent to solutions of the 
Bogoliubov-de Gennes equations [3]? Probably not. Al- 
though the ansatz in Eq. (4) is quite general and not re- 
stricted to spherical or axial symmetric condensates it al- 
lows for the description of the subset of excitations which 
are consistent with the ansatz (4). To obtain all solutions 
of the Bogoliubov-de Gennes equations with a variational 
approach the ansatz must be further generalized. This 
will be the objective of future work. 



VI. CONCLUSION 

We used an ansatz with several coupled Gaussian func- 
tions to obtain an improved variational description of 
the dynamics of Bose-Einstein condensates. We ap- 
plied a time-dependent variational principle to the Gross- 
Pitaevskii equation and obtained dynamical equations for 
the variational parameters with the improved variational 
method of coupled Gaussian functions. We discussed 
methods for solving these equations and analyzing the 
stability using methods of nonlinear dynamics. When 
can we expect the proposed variational methods to be 
better than a grid method? This is hard to say in general. 
The variational ansatz with a finite number of Gaussians 
is still an approximation. However, grid calculations are 
only "exact" in the limit of a small step size and thus an 
infinitely large grid. 

For a Gross-Pitaevskii equation with long-range 
monopolar (1/r) and dipolar (l/r^) interaction, the im- 
proved variational results, the convergence, and compar- 
isons with numerical calculations will be subject of the 
subsequent paper [19]. It will be shown that a low num- 
ber of Gaussians is sufficient to obtain converged results 
and thus the number of parameters in the variational 
computations is significantly smaller compared to calcu- 
lations on grids, and that a wealth of new phenomena 
is obtained by using the ansatz with coupled Gaussian 
functions. 
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